#Replication R script for tables and figures for 'Public Opinion and Nuclear Use'

#  This script is organized as a nested list
#  In RStudio, hit Alt+O

# ENVIRONMENT --------------------
##  Loading packages and data ----------

library(logr)
log_open()

rm(list=ls())
your_filepath = "C:/Users/m_h_g/Dropbox/Nuclear Survey Project/replication file/replication_bgg"
setwd(your_filepath)

library(tidyverse)
library(estimatr)
library(texreg)
library(xtable)
library(ggpubr)
#  Note: the plyr package is required to calculate the
#  bootstrap results, but it isn't loaded due to conflicts
#  with dplyr.

format_num <- function(x, digits = 3) {
  x <- as.numeric(x)
  return(paste0(sprintf(paste0("%.", digits, "f"), x)))
}
invlogit = function(x) 1 / (1+exp(-1*x))

#  script options for block bootstrap of model-based estimates ---
#  if fresh = F, precalculated results are loaded
#  if fresh = T, results are calculated from scratch
#  Nsim controls the number of bootstrap replicates
fresh = F
Nsim = 1000

##  Style for figures and tables -------

theme_allPlots =   
  function(theSize = 9){
    theme_bw() +
      theme(legend.title = element_blank(),
            axis.ticks = element_blank(),
            axis.title = element_text(size = theSize, color = "gray5"),
            axis.text =  element_text(size = theSize, color = "gray5"),
            legend.text = element_text(size = theSize, color = "gray5"),
            strip.text =  element_text(size = theSize, color = "gray5"),
            strip.background = element_blank(),
            axis.title.x = element_text(margin = margin(6,0,0,0)),
            axis.title.y = element_text(margin = margin(0,6,0,0)),
            plot.margin = margin(1,1,1,1),
            panel.grid = element_line(size = .25, color = "gray80")) 
  }


makezerovec = function(x) vapply(x, function(x) rep(0, as.numeric(x)) %>% paste(collapse=""), "a")
trail_zero  = function(x, digits_desired=2){
  digits = nchar(gsub(".+\\.", "", x))
  if(max(digits)>digits_desired) x = round(x, digits_desired)   # round it
  x = as.character(x)
  x[!grepl("\\.", x)] = paste0(x[!grepl("\\.", x)], ".")        # add dropped decimal points
  digits = nchar(gsub(".+\\.", "", x))
  paste0(x, makezerovec(digits_desired - digits))               # add dropped zeros
}

################################ ---------------
# LOADING DATA ------------------------

# Loading Main Study data -------
d_exper = read_csv("BGG_mainstudy_surveyExperiment.csv")

d_choice = read_csv("BGG_mainstudy_choiceExperiment.csv")

# Loading Pilot Study data ---------

s1_exper = read_csv("BGG_pilotexper.csv")

s1_conjoint = read_csv("BGG_pilot_conjoint_diffs_clean.csv")

# Loading PSV and ABP Replication data --------

abp = read_csv("abp.csv") 

psv = read_csv("psv.csv")

################################ ---------------
# EDITING DATA ------------------------
#   Mean imputation for main study -----------------

d_exper_meanImputed = d_exper %>%
  mutate_at(
    vars(ends_with("_N")),
    function(x){x[is.na(x)]=mean(x, na.rm=T); x-mean(x, na.rm=T)}) %>%
  mutate(
    educ = ifelse(is.na(educ), "Some college, no degree", educ),
    education = ifelse(is.na(education), "Some college, no degree", education),
    pol_interest = ifelse(is.na(pol_interest), "Somewhat interested", pol_interest),
    pol_discuss = ifelse(is.na(pol_discuss), "Once per week or more", pol_discuss)
  )

d_choice_meanImputed = d_choice %>%
  mutate_at(
    vars(ends_with("_N")),
    function(x){x[is.na(x)]=mean(x, na.rm=T); x-mean(x, na.rm=T)})%>%
  mutate(
    educ = ifelse(is.na(educ), "Some college, no degree", educ),
    education = ifelse(is.na(education), "Some college, no degree", education),
    pol_interest = ifelse(is.na(pol_interest), "Somewhat interested", pol_interest),
    pol_discuss = ifelse(is.na(pol_discuss), "Once per week or more", pol_discuss)
  )

#   Conjoint-style data for pilot study and main study ------------

d_conjoint = 
  bind_rows(
    d_choice %>% (function(x){names(x)=gsub("s2", "otherStrike", names(x)); x}) %>% mutate(strikeNum = 1),
    d_choice %>% (function(x){names(x)=gsub("s1", "otherStrike", names(x)); x}) %>% mutate(strikeNum = 2, chosen=1-chosen) %>% (function(x){names(x)=gsub("s2", "s1", names(x)); x}) %>%
      mutate_at(vars(starts_with("adv_")), function(x)x*-1) %>%
      mutate_at(vars(starts_with("disadv_")), function(x)x*-1)
  ) %>%
  mutate(types = paste(s1_strike, "vs.", otherStrike_strike),
         s1_chance = relevel(factor(s1_chance), "90 percent"),
         s1_civ = factor(s1_civ, c(10, 100, 1000)),
         s1_mil = recode_factor(s1_mil, No="None", Low="Low", High="High"),
         s1_env = factor(s1_env, c("Minimal", "Moderate", "Severe")),
         s1_strike = factor(s1_strike),
         s1_ally = recode_factor(s1_ally, Most="Most support", Few="Few support"))


#  make variables match

d_conjoint_s1 = mutate(s1_conjoint, adv_env=env_adv, adv_mil=milCas_adv, adv_follow=followUp_adv, adv_chance=chance_adv, adv_civ=civCas_adv, disadv_strike=1-nuke_adv)

d_conjoint_s2_for_merge = d_conjoint %>% 
  filter(numDisadv == "fully randomized") %>%
  mutate(
    ovr_adv = (adv_mil + adv_ally + adv_civ + adv_env + adv_chance),
  )

#  Merge study 1 and study 2 conjoints

d_conjoint_merged = bind_rows(
  d_conjoint_s1 %>% mutate(study=1, Study="Study 1") %>% select(-education, -ethnicity, -gender, -hhi, -hispanic, -political_party), 
  d_conjoint_s2_for_merge %>% mutate(study=2, Study="Study 2")
) %>%
  mutate(
    disadv_ally = -1*adv_ally,
    disadv_civ = -1*adv_civ,
    disadv_env = -1*adv_env,
    disadv_nuke = coalesce(-1*adv_strike, -1*nuke_adv)
  )

#   Recoding Aronow, Baron, and Pinson data ------------

abp = abp %>%
  mutate(
    educ = dplyr::recode(round(as.numeric(educ)), `1`="Did not complete high school", `2`="High school graduate",
                         `3`="Some college, no degree", `4`="Associate's degree",
                         `5`="Bachelor's degree", #`6`="Bachelor's degree", 
                         `6`="Graduate or professional degree"),
    
    exper_choice01 = pref,
    s_num = as.numeric(condition == "Nuclear Advantage"),
    d_num = 0,
    dataset = "abp",
    weight = weights
  ) %>%
  filter(!is.na(s_num))

#  Trim weights

abp$weight[abp$weight < .25] = .25
abp$weight[abp$weight > 4] = 4

#   Recoding Press, Sagan, and Valentino data --------------

psv = psv %>%
  mutate(
    educ = dplyr::recode(educ, `1`="Did not complete high school", `2`="High school graduate",
                         `3`="Some college, no degree", `4`="Associate's degree",
                         `5`="Bachelor's degree", `6`="Graduate or professional degree"),
    exper_choice01 = psv.pref,
    s_num = as.numeric(condition == "Nuclear Advantage"),
    d_num = 0,
    dataset = "psv",
    weight = weight
  ) %>%
  filter(!is.na(s_num))

#     Vignette experiment data, starting with Study 1 -----------
s1_exper = s1_exper %>%
  mutate(
    educ = education,
    d_num = 1,
    s_num = z,
    exper_choice01 = prefer_nuke_01,
    weight = 1,
    pol_interest = "Somewhat interested",
    pol_discuss = "Once per week or more"
  )

#Combine study 1 and study 2 experiment results (will be useful in making Figure
#D.4 in the Appendix)
s1s2_exper = bind_rows(
  s1_exper %>% mutate(study=1), 
  d_exper_meanImputed %>% mutate(study=2)
)%>%
  mutate_at(
    vars(age, ends_with("_N")),
    function(x){x[is.na(x)]=mean(x, na.rm=T); x-mean(x, na.rm=T)}) %>%
  mutate(
    dataset = "original",
    age2 = age^2,
    age3 = age^3
  )

#Combining Study 1, study 2, ABP, and PSV results - also for Figure D.4
all_exper =
  bind_rows(
    s1s2_exper,
    abp[,which(names(abp) %in% names(s1s2_exper))] %>% mutate(study=0),
    psv[,which(names(psv) %in% names(s1s2_exper))] %>% mutate(study=0) %>% select(-gender, -race)
  ) %>%
  mutate(
    age2 = age^2,
    age3 = age^3
  )

#     Add both treatments to merged conjoint; will be used in Table C.2
#     to make sure conjoint results are robust to treatment assignment ----------

d_conjoint_merged =
  left_join(
    d_conjoint_merged,
    all_exper %>% select(id, study, s_num, d_num)
  )

#     Bootstrap ---------

#  Note: you must install the plyr package
#  but plyr's large number of conflicts with dplyr
#  makes it inadvisable to load both packages at once.

llply = plyr::llply
alply = plyr::alply

expand_id_boot  = function(samp, id) map(samp, function(x)which(id %in% x)) %>% unlist

make_bootsample = function(id, n, seed=0){
  set.seed(seed)
  ids = replicate(n, sample(unique(id), replace=T)) 
  alply(ids, 2, (function(x)expand_id_boot(x, id)))
}

block_boot = function(df, fn, id_var="id", n, groups_in_output=NULL, returnBoots=F, needed = NULL, quantileType=7, seed=0){
  
  #  if specified, limit the data.frame to the variable names in "needed" (must be a character vector)
  if(!is.null(needed)) df = df[,needed]
  
  #  point estimates
  out = fn(df)
  
  #  which columns of the output will you need to match on when you merge the bootstrap replicates?
  if(is.null(groups_in_output)){groups = names(out)[1:(ncol(out)-1)]
  } else if(groups_in_output=="df"){groups = group_vars(df)
  } else if(groups_in_output=="out"){groups = group_vars(out)
  } else groups = groups_in_output
  boot = out[,groups,drop=F]   # has the output columns you'll need, and nothing else.
  
  boot_samples = make_bootsample(as.data.frame(df)[,id_var], n, seed)
  
  #  take each set of row numbers, create the data frame, run fn() on it, and merge the results.
  #     (joining this way makes sure that if for some reason, one of your groups is missing in one of your bootstrap replicates, it'll show up NA)
  booted = llply(boot_samples, (function(rownums) left_join(boot, fn(df[rownums,]), by=groups))) %>% map_df(function(x)t(x[,length(x)]))     
  booted = apply(booted, 2, as.numeric)
  
  #  now you have a data.frame of bootstrapped estimates.
  
  #  bind point estimates to the summary stats 
  out = out %>% cbind(
    se    = apply(booted, 1, sd, na.rm=T),
    q2.5  = apply(booted, 1, quantile, .025, na.rm=T, type=quantileType),
    q97.5 = apply(booted, 1, quantile, .975, na.rm=T, type=quantileType),
    gt0   = apply(booted, 1, function(x) mean(x>0, na.rm=T)),
    lt0   = apply(booted, 1, function(x) mean(x<0, na.rm=T)),
    NAs   = apply(booted, 1, function(x) sum(is.na(x)))
  ) %>%
    mutate(Nsim = ncol(booted))
  
  if(returnBoots) out = cbind(out, booted)
  
  out
}

#################################### ------------------
# RESULTS: VIGNETTE EXPERIMENT --------------------
#   CALC: Model estimates, both studies --------------------

#  Functions for computing model estimates
fn_logit_s1 = function(x){
  x %>% 
    do(
      glm(chosen ~ disadv_nuke + adv_chance + adv_mil + adv_follow + disadv_civ + disadv_env, 
          family = binomial(), data = .) %>% .$coefficients %>%
        (function(x){data.frame(variable = names(x), estimate = x)})
    )
}

fn_logit_s2 = function(x){
  x %>% 
    do(
      glm(chosen ~ disadv_nuke + adv_chance + adv_mil + disadv_ally + disadv_civ + disadv_env, 
          family = binomial(), data = .) %>% .$coefficients %>%
        (function(x){data.frame(variable = names(x), estimate = x)})
    )
}

fn_logit_restrict = function(x){
  x %>% 
    do(
      glm(chosen ~ disadv_nuke + ovr_adv, 
          family = binomial(), data = .) %>% .$coefficients %>%
        (function(x){data.frame(variable = names(x), estimate = x)})
    )
}

d_conjoint_merged %>% filter(study==1) %>% do(fn_logit_s1(.))
d_conjoint_merged %>% filter(study==2)%>% do(fn_logit_s2(.))
d_choice %>% filter(nukeAdv != "fully randomized") %>% do(fn_logit_s2(.))

#  If script is set to bootstrap fresh results:
#  Calculate parameter estimates that will go in the panels of Table E.1

if(fresh){
  Fit_s1 = block_boot(
    d_conjoint_merged %>% filter(study==1), fn_logit_s1, "id", Nsim
  )
  tab_modelFit_s2 = block_boot(
    d_conjoint_merged %>% filter(study==2), fn_logit_s2, "id", Nsim
  )
  
  tab_modelFit_restrict_s1 = block_boot(
    d_conjoint_merged %>% filter(study==1), fn_logit_restrict, "id", Nsim
  )
  tab_modelFit_restrict_s2 = block_boot(
    d_conjoint_merged %>% filter(study==2), fn_logit_restrict, "id", Nsim
  )
  
  #Restricted choice experiments parameter estimates
  tab_modelFit_restrict = block_boot(
    d_choice %>% filter(nukeAdv != "fully randomized"), fn_logit_s2, "id", Nsim
  )
  
  write_csv(
    bind_rows(
      tab_modelFit_restrict_s1 %>% mutate(study = 1),
      tab_modelFit_restrict_s2 %>% mutate(study = 2)
    ),
    "parametric_results_restricted.csv"
  )
  write_csv(
    bind_rows(
      tab_modelFit_s1 %>% mutate(study = 1),
      tab_modelFit_s2 %>% mutate(study = 2)
    ),
    "parametric_results.csv"
  )
}

#  If script is set to load precalculated results,
#  load them

if(!fresh){
  tab_modelFit = read_csv("parametric_results.csv")
  tab_modelFit_s1 = tab_modelFit %>% filter(study==1)
  tab_modelFit_s2 = tab_modelFit %>% filter(study==2)
  
  tab_modelFit_restrict = read_csv("parametric_results_restricted.csv")
  tab_modelFit_restrict_s1 = tab_modelFit_restrict %>% filter(study==1)
  tab_modelFit_restrict_s2 = tab_modelFit_restrict %>% filter(study==2)
}

#   TABLE 1: Study 2 Vignette Results, Full and Attentive Samples --------------------------

#  Fit all regressions

fit_full <-     lm_robust(exper_choice01 ~ s_num*d_num, data = d_exper_meanImputed, weights = weight)
fit_full_cov <- lm_robust(exper_choice01 ~ s_num*d_num + 
                            foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + foreign_preventNuke_N +
                            foreign_preventTerror_N + foreign_promoteDem_N + foreign_promoteHR_N + 
                            military_destroyTerrorist_N + military_helpUN_N + military_stopGenocide_N +
                            military_oil_N + military_promoteDem_N + military_protectAllies_N + 
                            age + pol_interest + education + pol_discuss, 
                          data = d_exper_meanImputed, weights = weight) 

fit_attentive     <-  lm_robust(exper_choice01 ~ s_num*d_num, data = d_exper_meanImputed, weights = weight,
                                subset = time_treatment > 60 & time_treatment < 600)
fit_attentive_cov <- lm_robust(exper_choice01 ~ s_num*d_num +
                                 foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + 
                                 foreign_preventNuke_N + foreign_preventTerror_N + foreign_promoteDem_N + 
                                 foreign_promoteHR_N + military_destroyTerrorist_N + military_helpUN_N + 
                                 military_stopGenocide_N + military_oil_N + military_promoteDem_N + 
                                 military_protectAllies_N + age + pol_interest + education + pol_discuss,
                               data = d_exper_meanImputed, weights = weight,
                               subset = time_treatment > 60 & time_treatment < 600)

# Create table

texreg(l = list(fit_full, fit_full_cov, fit_attentive, fit_attentive_cov),
       custom.model.names = c("No Controls", "Controls", "No Controls", "Controls"),
       custom.header = list("Full Sample" = 1:2, "Attentive Sample" = 3:4),
       custom.coef.names = c("$\\alpha$~~Constant", "\\\\[-2.5ex] $\\beta_1$~~Better chance of success", "\\\\[-2.5ex] $\\beta_2$~~More destructive", "\\\\[-2.5ex] $\\beta_3$~~Better chance $\\times$ more destructive"),
       omit.coef = "foreign|military|age|pol_|education",
       stars = c(0.02, 0.1, 0.2),
       digits = 3,
       caption.above = T,
       include.ci = F,
       custom.note = "OLS with robust SEs. Stars denote pre-registered one-tailed p-values below 0.01, 0.05, and 0.1.",
       caption = "Regression analysis of vignette experiment.",
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       label = "tab:attentive",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .) %>%
  write("C:/Users/m_h_g/Dropbox/Nuclear Survey Project/paper 2/tables/attentive_reg.tex")

#   Figure D.1: Study 2 Group Means -----------------------------------------------------

full <- d_exper_meanImputed %>%
  group_by(s_num, d_num) %>%
  do(
    lm_robust(exper_choice01 ~ 1, data = ., weights = weight) %>% tidy) %>%
  mutate(
    sample = "Full"
  )

attentive <- d_exper_meanImputed %>%
  group_by(s_num, d_num) %>%
  do(
    lm_robust(exper_choice01 ~ 1, data = ., weights = weight,
              subset = time_treatment >= 60 & time_treatment <= 600) %>% tidy) %>%
  mutate(
    sample = "Attentive"
  )


tab_plot <- rbind(full, attentive)
tab_plot$sample_f = factor(tab_plot$sample, levels=c('Full','Attentive'))
tab_plot <- tab_plot %>%
  mutate(
    effectiveness = ifelse(s_num == 0, "+0%", "+20%"),
    destructiveness = ifelse(d_num == 0, "Low Harm", "High Harm")
  )


tab_plot %>% ggplot(aes(x = effectiveness, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0) + 
  geom_path(aes(group = d_num, linetype = destructiveness)) +
  facet_wrap(~sample_f) + 
  ylim(-0.02, 0.55) +
  xlab("Nuclear military effectiveness (chance of success)") + 
  ylab("Prop. preferring nuclear") + 
  theme_bw() +
  theme(legend.title = element_blank(),
        legend.position = "bottom")# + 

ggsave("s2_attentive.pdf", height = 5, width = 7)
ggsave("FigD.1.tiff", height = 5, width = 7)


#   Figure D.2: S2 Robustness Charts ----------------------------------------------------

coef    <- rep(NA, 120)
se      <- rep(NA, 120)
pval    <- rep(NA, 120)
ci_low  <- rep(NA, 120)
ci_high <- rep(NA, 120)
cutoff  <- rep(NA, 120)

for(i in 45:120){
  fit        <- lm_robust(exper_choice01 ~ s_num*d_num, data = d_exper, 
                          weights = weight,
                          alpha = 0.1,
                          subset = time_treatment >= i & time_treatment <= 600)
  coef[i]    <- fit$coefficients[4]
  se[i]      <- fit$std.error[4]
  pval[i]    <- fit$p.value[4]
  ci_low[i]  <- fit$conf.low[4]
  ci_high[i] <- fit$conf.high[4]
  cutoff[i]  <- i
}

dat <- as_tibble(cbind(coef, se, pval, ci_low, ci_high, cutoff)) %>%
  slice(45:120)

ribbons <- dat %>% ggplot(aes(x = cutoff, y = coef)) +
  geom_path() +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "grey", alpha = 0.5) +
  geom_vline(xintercept = 60, color = "steelblue", linetype = "dashed") + 
  xlab("Minimum Seconds Spent to be Included") +
  ylab("Dif-in-Dif Coefficient (No Controls)") +
  ggtitle("(a) Maximum Seconds to be Included = 600") +
  theme_bw()

coef1    <- rep(NA, 120)
se1      <- rep(NA, 120)
pval1    <- rep(NA, 120)
ci_low1  <- rep(NA, 120)
ci_high1 <- rep(NA, 120)
cutoff1  <- rep(NA, 120)

for(i in 45:120){
  fit        <- lm_robust(exper_choice01 ~ s_num*d_num, data = d_exper, 
                          weights = weight,
                          alpha = 0.1,
                          subset = time_treatment >= i & time_treatment <= 480)
  coef1[i]    <- fit$coefficients[4]
  se1[i]      <- fit$std.error[4]
  pval1[i]    <- fit$p.value[4]
  ci_low1[i]  <- fit$conf.low[4]
  ci_high1[i] <- fit$conf.high[4]
  cutoff1[i]  <- i
}

dat1 <- as_tibble(cbind(coef1, se1, pval1, ci_low1, ci_high1, cutoff1)) %>%
  slice(45:120)

ribbons1 <- dat1 %>% ggplot(aes(x = cutoff1, y = coef1)) +
  geom_path() +
  geom_ribbon(aes(ymin = ci_low1, ymax = ci_high1), fill = "grey", alpha = 0.5) +
  geom_vline(xintercept = 60, color = "steelblue", linetype = "dashed") + 
  xlab("Minimum Seconds Spent to be Included") +
  ylab("Dif-in-Dif Coefficient (No Controls)") +
  ggtitle("(b) Maximum Seconds to be Included = 480") +
  theme_bw()

ribbon_panel <- gridExtra::arrangeGrob(ribbons, ribbons1, nrow = 2)

ggsave("ribbon_panel.pdf", ribbon_panel, height = 7, width = 5)
ggsave("FigD.2.tiff", ribbon_panel, height = 7, width = 5)

#   Figure D.3: Study 1 Group Means  --------------
s1full <- s1_exper %>%
  group_by(Z) %>%
  do(
    lm_robust(prefer_nuke_01 ~ 1, data = .) %>% tidy) %>%
  mutate(
    sample = "Full"
  )

s1attentive <- s1_exper %>%
  group_by(Z) %>%
  do(
    lm_robust(prefer_nuke_01 ~ 1, data = .,
              subset = time_treatment >= 60 & time_treatment <= 600) %>% tidy) %>%
  mutate(
    sample = "Attentive"
  )


s1tab_plot <- rbind(s1full, s1attentive)
s1tab_plot$sample_f = factor(s1tab_plot$sample, levels=c('Full','Attentive'))
s1tab_plot <- s1tab_plot %>%
  mutate(
    effectiveness = ifelse(Z == "Equally\neffective", "+0%", "+20%")
  )

s1tab_plot %>% ggplot(aes(x = effectiveness, y = estimate*100)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low*100, ymax = conf.high*100), width = 0) + 
  geom_path(aes(group = sample)) +
  ylim(0, 60) +
  facet_wrap(~sample_f) + 
  xlab("Nuclear military effectiveness (chance of success)") + 
  ylab("Percent preferring nuclear") + 
  theme_bw() 

ggsave("s1_attentive.pdf", height = 5, width = 7)
ggsave("FigD.3.tiff", height = 5, width = 7)

#All Vignette Experiment Results: Figure D.4 meta-analysis
#   FIGURE D.4: Meta-analysis -----------

tab_plot = 
  rbind(
    all_exper, #%>% filter(study!=0),
    all_exper %>% mutate(dataset = "Pooled", study =3)
  ) %>%
  filter(!is.na(s_num), !is.na(d_num), !is.na(dataset)) %>%
  group_by(s_num, d_num, dataset, study) %>%
  do(lm_robust(exper_choice01 ~ 1, data = ., weights = weight) %>% tidy)

tab_plot = tab_plot %>%
  group_by() %>%
  mutate(
    Dataset = recode_factor(dataset, abp = "Aronow et al.", psv="Press et al.", original="Original", Pooled="Pooled"),
    Facet = recode_factor(study, `0`="Extant studies", `1`="Pilot", `2`="Main Study", `3`="Meta-analysis")
  )

# ADD --- precision-weighted

tab_plot = tab_plot %>%
  mutate(Destruction = recode_factor(d_num, `1`="High Harm", `0`="Low Harm"),
         Effectiveness = recode_factor(s_num, `0`="+0%", `1`="+20%")) %>%
  mutate(Portrayal = Destruction)

g = tab_plot %>%
  ggplot(aes(x = Effectiveness, y = estimate, ymin = conf.low, ymax = conf.high, shape = Dataset, group = paste(Portrayal, Dataset))) +
  geom_errorbar(width = 0, lty = 1, position= position_dodge(width=.15), color = "gray5", size = .25) +
  geom_path(aes(lty = Portrayal), color = "gray50", position = position_dodge(width = .15)) +
  
  geom_point(size = 1.75, position= position_dodge(width=.15)) +
  #geom_rect(aes(xmin=NA_real_, ymin = NA_real_, xmax=NA_real_, ymax=NA_real_), color = "gray5", fill="white", size = .4) +
  
  theme_allPlots() +
  theme(
    legend.title = element_text(size = 9, hjust=0),
    axis.title.y = element_text(angle=0, vjust = .5),
    legend.margin = margin(-2,0,-2,0),
    panel.grid.major.x=element_blank()) +
  scale_shape_manual(values = c(24,22,21,#22,
                                23)) +
  scale_fill_manual(values=c("gray55", "gray95")) +
  scale_x_discrete(expand = c(.225, .225)) +
  scale_linetype_manual(values=c("solid", "dashed")) +
  scale_y_continuous(labels = (0:7)*10, breaks = (0:7)/10) +
  coord_cartesian(ylim = c(0, .65)) +
  guides(fill = guide_legend(rev=T)) +
  
  facet_wrap(~Facet, nrow=1) +
  labs(y = c("Percent\npreferring\nnuclear"),
       x = c("Nuclear military advantage (chance of success)"))

ggsave("surveyExper_vs_others.pdf", g, "pdf", width = 6.5, height = 2.75)   
ggsave("FigD.4.tiff", g, "tiff", width = 6.5, height = 2.75)

#   Table D.1: Regression meta-analysis ---------------------------

regs_all = list(
  lm_robust(exper_choice01 ~ s_num, 
            data = s1s2_exper %>% filter(study==1), weights = weight),
  lm_robust(exper_choice01 ~ s_num+ age + age2 + age3 +  
              educ, 
            data = s1s2_exper %>% filter(study==1), weights = weight),
  lm_robust(exper_choice01 ~ s_num*d_num, 
            data = s1s2_exper %>% filter(study==2), weights = weight),
  lm_robust(exper_choice01 ~ s_num*d_num+ 
              foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + foreign_preventNuke_N +
              foreign_preventTerror_N + foreign_promoteDem_N + foreign_promoteHR_N + 
              military_destroyTerrorist_N + military_helpUN_N + military_stopGenocide_N +
              military_oil_N + military_promoteDem_N + military_protectAllies_N + 
              age + age2 + age3 +pol_interest + educ + pol_discuss, 
            data = s1s2_exper %>% filter(study==2), weights = weight)
)

regs_meta = list(
  lm_robust(exper_choice01 ~ s_num*d_num, 
            data = s1s2_exper, weights = weight),
  lm_robust(exper_choice01 ~ s_num*d_num+  age + age2 + age3 +
              educ , 
            data = s1s2_exper, weights = weight),
  lm_robust(exper_choice01 ~ s_num*d_num, data = all_exper),
  lm_robust(exper_choice01 ~ s_num*d_num + age + age2 + age3 +educ, data = all_exper)
)
out = texreg(c(regs_all, regs_meta),
             #custom.model.names = c("No controls", "Controls"),
             custom.coef.map = list(`(Intercept)` = "$\\alpha$~~Constant", 
                                    s_num = "$\\beta_1$~~Better chance of success", 
                                    d_num = "$\\beta_2$~~More destructive", 
                                    `s_num:d_num` = "$\\beta_3$~~Better chance $\\times$ more destructive"),
             omit.coef = "foreign|military|age|pol_|educ",
             # reorder.coef = c(3,2,4,1),
             booktabs = T,
             caption.above = T,
             stars = c(0.02, 0.1, 0.2),
             digits = 3,
             custom.note = "OLS estimates. Robust SEs in parentheses. Controls and one-tailed tests pre-registered. $^*p<0.1$, $^{**}p<0.05$, $^{***}p<0.01$.",
             # star.symbol = c("*", "**"),
             caption = "Conditional effects of strike attributes.",
             include.rsq = F,
             include.ci=F,
             include.rmse =F,
             use.packages = F) 

out = strsplit(out, "\n")[[1]]
out[11] = gsub("of success", "",         out[11])
out[12] = paste("~~~~~of success",       out[12])
out[15] = gsub("more destructive", "",   out[15])
out[16] = paste("~~~~~more destructive", out[16])

out = ifelse(grepl("beta", out), paste("\\\\[-1.5ex]", out), out)

out = paste(out[9:length(out)], collapse = "\n")

out %>% gsub(".+.Model 8", "", .) %>%
  gsub(".end.tabular.+", "", .) %>%
  gsub("p<0.1", "p<0.05", .) %>%
  gsub("p<0.02", "p<0.01", .) %>%
  write("C:/Users/m_h_g/Dropbox/Nuclear Survey Project/paper 2/tables/tab_regs_pooled.txt")


#################################### -----------------
# RESULTS: CONJOINT ------------
#   Figure 1: Illustrative example --------

#  Chance of success facet
tab_chance <- d_choice %>%
  filter(Randomization=="Realistic") %>%
  filter(s1_strike == "Nuclear", s2_strike == "Conventional") %>%
  filter(disadv_index == 3 | disadv_index == 0) %>%
  filter(nukeAdv == "neither" | nukeAdv == "success") %>%
  mutate(side_effects = case_when(disadv_index == 3 ~ "High Harm", T ~ "Low Harm"),
         chance = case_when(nukeAdv == "success" ~ "+20% chance", T ~ "+0% chance")) %>%
  group_by(chance, side_effects, disadv_index, nukeAdv) %>%
  do(
    lm_robust(chosen ~ 1, data = ., clusters = id, weights = weight) %>% tidy
  )
tab_chance$plot <- rep("Chance", 4)

chance_plot <- tab_chance %>% ggplot(aes(x = reorder(chance, estimate), y = estimate*100)) +
  geom_point() +
  geom_line(aes(group = side_effects, lty = side_effects)) +
  geom_errorbar(aes(ymin=conf.low*100, ymax=conf.high*100), width=0, position = position_dodge(width=0)) +
  ylim(10, 70) + 
  theme_allPlots() +
  ylab("Percent \npreferring \nnuclear") +
  xlab("Nuclear military advantage") + 
  ggtitle("Effect of chance of success") +
  theme(legend.position = "none",
        axis.title.y = element_text(angle = 0, vjust=.5, hjust=.5),
        plot.title = element_text(size = 9, hjust = 0.5))

#     Casualties facet ----------

tab_mil_cas2 <- d_choice %>%
  filter(Randomization=="Realistic") %>%
  filter(s1_strike == "Nuclear", s2_strike == "Conventional") %>%
  filter(disadv_index == 3 | disadv_index == 0) %>%
  filter(nukeAdv == "neither" | nukeAdv == "military") %>%
  mutate(side_effects = case_when(disadv_index == 3 ~ "High Harm", T ~ "Low Harm"),
         mil_cas = case_when(nukeAdv == "military" ~ "Fewer casualties", T ~ "Same casualties") %>% factor(unique(.)[2:1])) %>%
  group_by(mil_cas, side_effects, disadv_index, nukeAdv) %>%
  do(
    lm_robust(chosen ~ 1, data = ., clusters = id, weights = weight) %>% tidy
  )

tab_mil_cas2$plot <- rep("U.S. Military Casualties", 4)

mil_cas_plot2 <- tab_mil_cas2 %>% ggplot(aes(x = reorder(mil_cas, estimate), y = estimate*100)) +
  geom_point() +
  geom_line(aes(group = side_effects, lty = side_effects)) +
  geom_errorbar(aes(ymin=conf.low*100, ymax=conf.high*100), width=0, position = position_dodge(width=0)) +
  ylim(10, 70) + 
  theme_allPlots() +
  ylab("") +
  xlab("Nuclear military advantage") + 
  ggtitle("Effect of military casualties") + 
  theme(legend.position = "none",
        plot.title = element_text(size = 9, hjust = 0.5))

mil_cas_plot2 

#     Combine and export, produces Figure 1
side_by_side_two2 <- ggarrange(chance_plot, mil_cas_plot2, nrow = 1, common.legend = TRUE, legend="bottom", widths = c(.45, .4))
ggsave("side_by_side_two2.pdf", side_by_side_two2, width = 6, height = 4)
ggsave("Fig1.tiff", side_by_side_two2, width = 6, height = 4)

#     Creating Appendix Table C.1, Means for Figure 1

out = bind_rows(
  tab_chance %>% dplyr::rename(Scenario=chance) %>% mutate(panel=1),
  tab_mil_cas2 %>% dplyr::rename(Scenario=mil_cas) %>% mutate(panel=2)
) %>%  
  group_by() %>%
  mutate(
    CI   = paste0("(", format_num(conf.low), ", ", format_num(conf.high), ")")
  ) %>%
  arrange(panel, side_effects, Scenario) %>%
  select(
    `Side effects`=side_effects, Scenario, Estimate=estimate, SE=std.error, CI
  )

for(i in nrow(out):2){
  if(out$`Side effects`[i]==out$`Side effects`[i-1]) out$`Side effects`[i] = ""
}

out %>%
  xtable(digits=3) %>%
  print(hline.after = nrow(.), only.contents=T, include.rownames=F, include.colnames = F) %>%
  write("tab_fig3_numeric.txt")

#   Restricted Choice Experiment Results, Table 2 -------------------------


lm_naked = d_choice %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv + adv_chance*numDisadv, data = ., clusters = id, weights = weight)# %>% tidy %>% mutate_at(vars(-term, -outcome), round, 3)
lm_covar = d_choice_meanImputed %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv + adv_chance*numDisadv  + 
              foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + foreign_preventNuke_N +
              foreign_preventTerror_N + foreign_promoteDem_N + foreign_promoteHR_N + 
              military_destroyTerrorist_N + military_helpUN_N + military_stopGenocide_N +
              military_oil_N + military_promoteDem_N + military_protectAllies_N + age + pol_interest + education, 
            data = ., clusters = id
  )

library(texreg)

texreg(list(lm_naked, lm_covar), 
       digits = 3,
       custom.coef.names = c("$\\alpha_0$~~Constant", "$\\delta_0$~~Disadvantages (0-3 scale)", "$\\alpha_M$~~Fewer military casualties", "$\\alpha_S$~~Better chance of success",
                             "$\\delta_M$~~Disadvantages $\\times$ fewer mil. casualties", "$\\delta_S$~~Disadvantages $\\times$ better chance"),
       custom.model.names = c("No controls", "Controls"),
       omit.coef = "foreign|military|age|interest|education",
       booktabs = T,
       caption.above = T,
       stars = c(0.02, 0.1, 0.2),#numeric(0),
       caption = "Conditional effects of strike attributes.",
       custom.note = "",
       include.rsq = F,
       include.ci = F,
       include.rmse =F) %>% 
  gsub(".+begin.center.|.label.+", "", .) %>%
  write("C:/Users/m_h_g/Dropbox/Nuclear Survey Project/paper 2/figures/tab_restricted_choice.txt")


#   APPENDIX TABLE C.2:  Controlling for survey experiment treatment status ---------------------

lm_naked = d_choice %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv + adv_chance*numDisadv + s_num*d_num, data = ., clusters = id, weights = weight)# %>% tidy %>% mutate_at(vars(-term, -outcome), round, 3)

lm_covar = d_choice_meanImputed %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv + adv_chance*numDisadv  + 
              foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + foreign_preventNuke_N +
              foreign_preventTerror_N + foreign_promoteDem_N + foreign_promoteHR_N + 
              military_destroyTerrorist_N + military_helpUN_N + military_stopGenocide_N +
              military_oil_N + military_promoteDem_N + military_protectAllies_N + age + pol_interest + education+ s_num*d_num, 
            data = ., clusters = id
  )

out = texreg(list(lm_naked, lm_covar), 
             digits = 3,
             custom.coef.map = list(`(Intercept)` = "$\\alpha_0$~~Constant", numDisadv = "$\\delta_0$~~Disadvantages (0-3 scale)", 
                                    adv_mil = "$\\alpha_M$~~Fewer military casualties", adv_chance = "$\\alpha_S$~~Better chance of success",
                                    `numDisadv:adv_mil` = "$\\delta_M$~~Disadvantages $\\times$ fewer mil. casualties", `numDisadv:adv_chance` = "$\\delta_S$~~Disadvantages $\\times$ better chance",
                                    s_num = "Better chance of success", d_num = "Equal destruction", `s_num:d_num` = "Better chance $\\times$ equal destruct."),
             custom.model.names = c("No controls", "Controls"),
             omit.coef = "foreign|military|age|interest|education",
             booktabs = T,
             caption.above = T,
             stars = numeric(0),#c(0.01, 0.05),
             caption = "Conditional effects of strike attributes.",
             custom.note = "",
             include.rsq = F,
             include.ci = F,
             include.rmse =F) %>% 
  gsub(".+begin.center.|.label.+", "", .) %>%
  strsplit("\n") %>%
  .[[1]]

for(i in 1:length(out)){
  if(grepl("alpha|delta|Better|Equal", out[i]) & !grepl("Constant", out[i])) out[i-1] = paste(out[i-1], " \\\\[-2ex]")
}

write(out, "C:/Users/m_h_g/Dropbox/Nuclear Survey Project/paper 2/figures/tab_restricted_choice_appendix.txt")

#   APPENDIX TABLE C.3:  Interacting survey experiment treatment status w/ parameters of interest -------------------

lm_naked = d_choice %>%
  mutate(Z=paste0(s_num, d_num)) %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv*Z + adv_chance*numDisadv*Z, data = ., clusters = id, weights = weight)# %>% tidy %>% mutate_at(vars(-term, -outcome), round, 3)

lm_covar = d_choice_meanImputed %>%
  mutate(Z=paste0(s_num, d_num)) %>%
  filter(nukeAdv != "fully randomized") %>%
  mutate(numDisadv = as.numeric(numDisadv)) %>%
  lm_robust(chosen ~ numDisadv + adv_mil + adv_chance + adv_mil*numDisadv*Z + adv_chance*numDisadv*Z +
              foreign_defendAllies_N + foreign_energySupply_N + foreign_helpUN_N + foreign_preventNuke_N +
              foreign_preventTerror_N + foreign_promoteDem_N + foreign_promoteHR_N + 
              military_destroyTerrorist_N + military_helpUN_N + military_stopGenocide_N +
              military_oil_N + military_promoteDem_N + military_protectAllies_N + age + pol_interest + education, 
            data = ., clusters = id
  )

out = texreg(list(lm_naked, lm_covar), 
             digits = 3,
             custom.coef.map = list(`(Intercept)` = "$\\alpha_0$~~Constant", `Z01`=" . . . $\\times$ ($d=0, s=1$)", `Z10`=" . . . $\\times$ ($d=1, s=0$)", `Z11`=" . . . $\\times$ ($d=1, s=1$)",
                                    numDisadv = "$\\delta_0$~~Disadvantages (0-3 scale)", `numDisadv:Z01`=" . . . $\\times$ ($d=0, s=1$) ", `numDisadv:Z10`=" . . . $\\times$ ($d=1, s=0$) ", `numDisadv:Z11`=" . . . $\\times$ ($d=1, s=1$) ",
                                    adv_mil = "$\\alpha_M$~~Fewer military casualties", `adv_mil:Z01`=" . . . $\\times$ ($d=0, s=1$)  ", `adv_mil:Z10`=" . . . $\\times$ ($d=1, s=0$)  ", `adv_mil:Z11`=" . . . $\\times$ ($d=1, s=1$)  ",
                                    adv_chance = "$\\alpha_S$~~Better chance of success", `adv_chance:Z01`=" . . . $\\times$ ($d=0, s=1$)   ", `adv_chance:Z10`=" . . . $\\times$ ($d=1, s=0$)   ", `adv_chance:Z11`=" . . . $\\times$ ($d=1, s=1$)   ",
                                    `numDisadv:adv_mil` = "$\\delta_M$~~Disadvantages $\\times$ fewer mil. casualties", `numDisadv:adv_mil:Z01`=" . . . $\\times$ ($d=0,s=1$)", `numDisadv:adv_mil:Z10`=" . . . $\\times$ ($d=1,s=0$)", `numDisadv:adv_mil:Z11`=" . . . $\\times$ ($d=1,s=1$)",
                                    `numDisadv:adv_chance` = "$\\delta_S$~~Disadvantages $\\times$ better chance", `numDisadv:adv_chance:Z01`=" . . . $\\times$ ($d=0,s=1$) ", `numDisadv:adv_chance:Z10`=" . . . $\\times$ ($d=1,s=0$) ", `numDisadv:adv_chance:Z11`=" . . . $\\times$ ($d=1,s=1$) "#,
                                    #s_num = "Better chance of success", d_num = "Equal destruction", `s_num:d_num` = "Better chance $\\times$ equal destruct."
             ),
             custom.model.names = c("No controls", "Controls"),
             omit.coef = "foreign|military|age|interest|education",
             booktabs = T,
             caption.above = T,
             stars = numeric(0),#c(0.01, 0.05),
             caption = "Conditional effects of strike attributes.",
             custom.note = "",
             include.rsq = F,
             include.ci = F,
             include.rmse =F) %>% 
  gsub(".+begin.center.|.label.+", "", .) %>%  
  strsplit("\n") %>%
  .[[1]]

for(i in 1:length(out)){
  if(grepl("alpha|delta", out[i]) & !grepl("Constant", out[i])) out[i-1] = paste(out[i-1], " \\\\[-1.5ex]")
}

write(out, "C:/Users/m_h_g/Dropbox/Nuclear Survey Project/paper 2/figures/tab_restricted_choice_appendix_interact.txt")

#   FIGURE D.5: Restricted choice means --------------------

tab_restrictedMeans = d_choice %>%
  filter(nukeAdv != "fully randomized") %>%
  filter(!is.na(nukeAdv), !is.na(numDisadv)) %>%
  filter(adv_chance>=0, adv_mil>=0) %>%
  group_by(nukeAdv, numDisadv) %>% 
  do(
    lm_robust(chosen ~ 1, data = ., clusters = id, weights = weight) %>% tidy
  )

tab_restrictedMeans = tab_restrictedMeans %>%
  group_by() %>%
  mutate(
    Label = recode_factor(nukeAdv, both = "Fewer military casualties,\nBetter chance of success", military = "Fewer military casualties,\nEqual chance of success",
                          success = "Equal military casualties,\nBetter chance of success", neither = "Equal military casualties,\nEqual chance of success") %>% factor(levels(.)[4:1]),
    Label2 = recode_factor(nukeAdv, 
                           both = "Fewer casualties, better chance", 
                           military = "Fewer casualties, equal chance",
                           success = "Equal casualties, better chance",#\n(most similar to BGG s=1)", 
                           neither = "Equal casualties, equal chance") %>% #\n(most similar to BGG s=0)") %>% 
      factor(levels(.)[4:1]),
    numDisadv = as.numeric(numDisadv)
  )

tab_plot = 
  tab_restrictedMeans %>% 
  mutate(Estimates = "Non-parametric") 

rightBorder= 0.27

g = tab_plot %>%
  ggplot(aes(x = -1*numDisadv, y = estimate, color = Estimates, shape = Label, fill = Label,
             lty = Estimates, size = Estimates,
             ymin = conf.low, ymax = conf.high, group = grepl("Equal ch", Label))) +
  geom_hline(aes(yintercept = .5), color = "white", size = .25) +
  geom_hline(aes(yintercept = .5), color = "gray75", size = 1.25, lty=5) +
  
  geom_rect(aes(xmin = rightBorder, xmax = Inf, ymin = -Inf, ymax = Inf), color = "white", fill = "white", lty=1, size=.25) +
  geom_rect(aes(xmin = -Inf, xmax = rightBorder, ymin = -Inf, ymax = Inf), color = "gray5", fill = "transparent", lty=1, size=.25) +
  
  # geom_rect(aes(xmin = .21, xmax = -.25, ymin = .13, ymax = .54), fill = "gray75", color = "gray55", size = .25, show.legend = F, lty = 1, alpha=.02) +
  # geom_text(aes(x = -.22, y = .338, label = "Most\nsimilar\nto PSV,\nABP, and\nBGG d=0"), color = "gray5", size = 2, hjust=0) +
  # 
  # geom_rect(aes(xmin = -3.18, xmax = -2.25, ymin = .115, ymax = .29), fill = "gray75", color = "gray55", size = .25, show.legend = F, lty = 1, alpha=.02) +
  # geom_text(aes(x = -2.85, y = .201, label = "Most similar\nto BGG d=1"), color = "gray5", size = 2, hjust=0) +
  
  
  geom_text(aes(x = 0.35, label = Label2, y = ifelse(numDisadv==0, estimate, NA_real_)), hjust = 0, size = 3.2,color = "gray5") +
  geom_path(aes(group = paste(Label, Estimates)), position = position_dodge(width=.15)) +
  geom_errorbar(width = 0, lty = 1, position= position_dodge(width=.15), color = "gray5", size = .4) +
  geom_point(
    aes(y = ifelse(Estimates=="Non-parametric", estimate, NA_real_)),
    size = 2.5, position= position_dodge(width=.15), color = "gray5",
    show.legend = F
  ) +
  scale_shape_manual(values=c(22,24,22,24)) +    #  FILL, THESE IN
  scale_fill_manual(values=c("gray60", "gray60", "gray95", "gray95")) +
  scale_size_manual(values = c(.5, 2)) +
  scale_linetype_manual(values = c(2,2)) +
  scale_color_manual(values = c("gray10", "gray85")) +
  scale_x_continuous(breaks = -3:0, labels = 3:0) +
  scale_y_continuous(breaks = (1:7)/10, labels = (1:7)*10) +
  coord_cartesian(xlim = c(-3, 2.15), ylim = c(.1, .7)) +
  guides(shape = guide_legend(include = F)) +
  theme_allPlots() +
  theme(#legend.position = c(.16, .93),#
    legend.position = "none",
    legend.key.width = unit(.6, "cm"),
    legend.key.height = unit(.45, "cm"),
    legend.background = element_blank(),
    #legend.title = element_text(size = 9, hjust=0),
    panel.border = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    
    axis.title.x = element_text(hjust = .21),
    axis.title.y = element_text(angle=0, vjust = .5, margin = margin(0,10,0,0)),
    panel.grid.minor.x = element_blank()) +
  labs(x = "Number of nuclear disadvantages", y = "Percent\npreferring\nnuclear")

ggsave("restricted_means.pdf", g, "pdf", width = 5, height = 3, dpi=1000)
ggsave("FigD.5.tiff", g, "tiff", width = 5, height = 3, dpi=1000)

#       Appendix Table D.2, showing estimates in Figure D.5 w/ numbers -----------------

out = tab_restrictedMeans %>% 
  mutate(CI = paste0("(", round(conf.low, 2), ", ", round(conf.high, 2), ")"),
         numDisadv = as.character(numDisadv),
         Label = gsub("\\n", " ", Label)) %>%
  select(
    Advantages = Label,
    Disadvantages = numDisadv,
    Estimate = estimate, `S.E.`=std.error, CI) %>%
  arrange(Advantages, desc(Disadvantages))

for(i in nrow(out):2){
  if(out$Advantages[i] == out$Advantages[i-1]) out$Advantages[i] = ""
}

for(i in nrow(out):2){
  if(out$Advantages[i-1] != ""){
    out$Advantages[i] = gsub(".+, " , "", out$Advantages[i-1])
    out$Advantages[i-1] = gsub(", .+" , ",", out$Advantages[i-1])
  }
}

out %>%
  xtable(caption = "Group means, restricted choice experiment, Study 2", label = "tab:restrictedMeans") %>%
  print(include.rownames = F, caption.placement = "top", booktabs = T, size = "\\small") %>%
  write("study2_restricted_means.txt")

#   Figure D.7: Main Study AMCE-------------------------

tab_amce = d_conjoint %>%
  filter(numDisadv == "fully randomized") %>%
  lm_robust(
    chosen ~ s1_strike + s1_chance + s1_mil + s1_ally + s1_civ + s1_env,
    data = ., clusters = id, weights = weight
  ) %>%
  tidy %>%
  filter(term!="(Intercept)")

tab_amce = tab_amce %>%
  mutate(value = gsub("s1_strike|s1_chance|s1_mil|s1_ally|s1_civ|s1_env", "", term))

# theLevels = levels(tab_amce$Variable)

tab_amce = 
  left_join(
    d_conjoint %>% select(s1_strike, s1_chance, s1_mil, s1_ally, s1_civ, s1_env) %>%
      gather(variable, value) %>% distinct,
    tab_amce
  )

tab_amce = tab_amce %>%
  group_by() %>%
  mutate(
    Variable = substr(variable,4,6),
    Variable = recode_factor(Variable, 
                             str = "Strike\ntype",
                             cha = "Chance of\nsuccess",
                             civ = "Civilian\ncasualties\n",
                             env = "Physical\ndestriction\n",
                             mil = "Military\ncasualties\n",
                             all = "Support of\nallies")
  )


tab_amce = tab_amce %>%
  mutate_at(vars(estimate, conf.low, conf.high), function(x)ifelse(!is.na(x), x, runif(length(x), -.000001, .000001)))

tab_amce = tab_amce %>% arrange(Variable, estimate) %>% mutate(value = factor(value, unique(value)))

g = tab_amce %>%
  ggplot(aes(y = estimate, x = value)) +
  
  geom_hline(aes(yintercept = 0), size = .5, color = "white") +
  geom_hline(aes(yintercept = 0), size = .5, color = "gray60", lty = 2) +
  
  geom_point(size = 1.5) +
  
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = -.31), fill = "white", color = "white") +
  geom_text(aes(y = -.409, label = value), size = 3.2, hjust = 0, color = "gray5") +
  
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = .5) +
  
  facet_wrap(~Variable, scales = "free_y", ncol = 1, strip.position = "left") +
  coord_flip() +
  theme_allPlots(9) +
  theme(strip.placement = "outside", 
        strip.text.y.left = element_text(angle = 0, hjust = 0, color = "gray5", margin=margin(0,3,0,0)),
        axis.text.y = element_blank(),
        #axis.title = element_blank(),
        axis.title.y = element_text(margin = margin(0,10,0,0)),
        panel.grid.minor.x = element_blank(),
        plot.margin = margin(2),
        panel.spacing = unit(.2, "cm"),
        panel.border = element_blank(),
        strip.background = element_blank()) +
  scale_y_continuous(expand = c(0, 0), limits = c(-0.41, 0.01), breaks = (-8:1)/20, labels = c("", "", (-6:1)/20)) +
  geom_hline(aes(yintercept = .7), color = "white", size = .5) +
  geom_hline(aes(yintercept = .75), color = "gray92", size = .25)+
  geom_hline(aes(yintercept = .85), color = "gray92", size = .25)+
  geom_hline(aes(yintercept = .95), color = "gray92", size = .25)+
  labs(x = "Attribute", y = "                 Effect on Pr(choose strike option)")


library(grid)
g1 = ggplotGrob(g)
g1$heights[as.character(g1$heights) == "1null"] = unit(c(2, 2, 3, 3, 3, 2), "null")
grid.newpage()
grid.draw(g1)

ggsave("amce.pdf", g1, "pdf", width = 5, height = 3.3)
ggsave("FigD.7.tiff", g1, "tiff", width = 5, height = 3.3)

#   Figure D.6: Pilot Study AMCE ---------------------------

d_amce = 
  s1_conjoint %>%
  select(-chosen_1) %>%
  (function(x){names(x) = gsub("_1", "", names(x)); x}) %>%
  mutate(
    milCas = factor(milCas, c("Minimal", "Moderate", "High")),
    env    = factor(env,    c("Minimal", "Moderate", "High")),
    civCas = factor(civCas),
    chance = factor(paste(chance, "percent"), c("90 percent", "70 percent")),
    followUp = factor(followUp, c("Yes", "No"))
  )

tab_amce = d_amce %>%
  lm_robust(chosen ~ followUp + chance + civCas + env + type + milCas, data = .) %>%
  tidy %>%
  mutate(
    first3   = substr(term, 1, 3),
    value = gsub("followUp|chance|civCas|env|type|milCas", "", term)
  )

tab_amce = 
  left_join(
    d_amce %>% 
      select(followUp, chance, civCas, env, type, milCas) %>%
      gather(variable, value) %>%
      distinct %>%
      mutate(first3 = substr(variable, 1, 3)),
    tab_amce
  ) %>%
  mutate_at(
    vars(estimate),
    function(x){x[is.na(x)] = runif(sum(is.na(x)), -.00001, .00001); x} 
  ) %>%
  mutate(
    category = recode(first3, 
                      fol = "Follow-up\nstrike", 
                      cha = "Chance of\nsuccess",
                      civ = "Civilian\ncasualties\n",
                      env = "Environmental\ndamage\n",
                      typ = "Strike\ntype",
                      mil = "Military\ncasualties\n")
  )

tab_amceMean = 
  d_amce %>% 
  select(followUp, chance, civCas, env, type, milCas, chosen) %>%
  gather(variable, value, -chosen) %>%
  group_by(variable, value) %>%
  filter(!is.na(value)) %>%
  summarize(mu = mean(chosen, na.rm=T),
            Percent = round(100*mu, 1))

g = left_join(
  tab_amce,
  tab_amceMean
) %>%
  mutate(
    value = ifelse(first3=="mil", paste0(value, " "), value),
    first3 = factor(first3, c("typ", "fol", "cha", "civ", "env", "mil"))
  ) %>%
  arrange(first3, estimate) %>%
  mutate(value = factor(value, value), percent=Percent/100) %>%
  ggplot(aes(y = estimate, x = value)) +
  
  geom_hline(aes(yintercept = 0), size = .5, color = "white") +
  geom_hline(aes(yintercept = 0), size = .5, color = "gray60", lty = 2) +
  
  geom_point(shape = 15) +
  
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = -.27), fill = "white", color = "white") +
  geom_text(aes(y = -.355, label = value), size = 3.5, hjust = 0, color = "gray5") +
  
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, size = .25) +
  
  facet_wrap(~category, scales = "free_y", ncol = 1, strip.position = "left") +
  coord_flip() +
  theme_minimal() +
  theme(strip.placement = "outside", 
        strip.text.y.left = element_text(angle = 0, hjust = 0, color = "gray5", size = 10),
        axis.text.y = element_blank(),
        axis.text.x = element_text(size = 9, color = "gray5"),
        axis.title = element_text(size = 10, color = "gray5"),
        axis.title.x = element_text(hjust = .7),
        #axis.title = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.spacing = unit(.35, "cm"),
        strip.background = element_blank()) +
  scale_y_continuous(expand = c(0, 0), limits = c(-0.36, 0.07), breaks = (-7:1)/20, labels = c("", "", (-5:1)/20)) +
  geom_hline(aes(yintercept = .7), color = "white", size = .5) +
  geom_hline(aes(yintercept = .75), color = "gray92", size = .25)+
  geom_hline(aes(yintercept = .85), color = "gray92", size = .25)+
  geom_hline(aes(yintercept = .95), color = "gray92", size = .25)+
  labs(x = "Attribute", y = "Effect on Pr(choose strike option)")

g

library(grid)
g1 = ggplotGrob(g)
g1$heights[as.character(g1$heights) == "1null"] = unit(c(2, 3, 3, 2, 3, 2), "null")
grid.newpage()
grid.draw(g1)


#out

ggsave("Study1_amceonly.png", g1, "png", width = 5.9, height = 3.6, dpi = 500)
ggsave("FigD.6.tiff", g1, "tiff", width = 5.9, height = 3.6, dpi = 500)

#Figure D.8: Net advantage figure
tab_plot_ovr_index = d_conjoint_merged %>%
  group_by(Study, ovr_adv, types) %>%
  do(lm_robust(chosen ~ 1, data = ., clusters = id) %>% tidy) %>%
  filter(types %in% c("Nuclear vs. Conventional", "Conventional vs. Conventional"))

tab_plot_ovr_index_N = d_conjoint_merged %>%
  group_by(Study, ovr_adv, types) %>%
  dplyr::summarize(N=n())


g = tab_plot_ovr_index %>%
  filter(types %in% c("Nuclear vs. Conventional", "Conventional vs. Conventional")) %>%
  group_by() %>%
  mutate(Study = case_when(Study == "Study 1" ~ "Pilot",
                           Study == "Study 2" ~ "Main Study")) %>%
  ggplot(aes(x = ovr_adv, y = estimate, ymin = conf.low, ymax = conf.high, color = types, lty = types, group=types)) +
  geom_point(shape=15) +
  geom_errorbar(width = .1, lty = 1, size = .25, alpha = .75) +
  geom_path() +
  
  scale_color_manual(values = c("gray15", "gray15")) +
  scale_linetype_manual(values = c(2, 1)) +
  scale_x_continuous(breaks = -5:5) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 1)) +
  geom_text(aes(x=-2, y=ifelse(ovr_adv==-4, .625, NA_real_), label="Conventional vs.\nConventional"), size=2.8) +
  geom_text(aes(x=2.5, y=ifelse(ovr_adv==-4, .375, NA_real_), label="Nuclear vs.\nConventional"), size=2.8) +
  
  theme_allPlots() +
  theme(
    legend.position = "none",
    legend.title = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.spacing.x = unit(.5, "cm")
  ) +
  labs(x = "Net number of disadvantages (negative) or advantages (positive)",
       y = "Percent choosing") +
  facet_wrap(~Study, scales="free_x")

ggsave("ovr_adv_both_conjoints.pdf", g, "pdf", width=6.5, height=3.5)
ggsave("FigD.8.tiff", g, "tiff", width=6.5, height=3.5)
##Discrete choice model parameter estimates: Figure E.1 and Table E.1
#   Appendix table E.1, unrestricted choice experiment parameter estimates ------------

#Making panel A of Table E.1 ("Pilot Study, Unrestricted)
tab_modelFit_s1 %>%
  filter(variable != "(Intercept)") %>%
  mutate(
    Variable = recode_factor(variable, disadv_nuke = "eta Nuclear", adv_chance = "alpha1 Chance of success", adv_mil = "alpha2 Prevents casualties", adv_follow = "alpha3 Follow-up strike",
                             disadv_civ = "delta1 Civilian casualties", disadv_env = "delta2 Environmental damage"),
    Estimate = paste0(format(round(estimate, 2)), " (", format(round(se, 2)), ")")
  ) %>%
  arrange(Variable) %>%
  select(Variable, Estimate) %>%
  xtable() %>%
  print(only.contents=T, hline.after = nrow(.),include.rownames=F, include.colnames=F) %>%
  gsub("eta", "$\\\\eta$~ & ", .)%>%
  gsub("alpha1", "$\\\\beta_1$~~ & ", .)%>%
  gsub("alpha2", "$\\\\beta_2$~~ & ", .)%>%
  gsub("alpha3", "$\\\\beta_3$~~ & ", .)%>%
  gsub("delta1", "$\\\\beta_4$~~ & ", .)%>%
  gsub("delta2", "$\\\\beta_5$~~ & ", .)%>%
  # gsub("delta3", "$\\\\beta_6$~~ & ", .)%>%
  write("tab_model_s1.txt")

#Making Panel B of Table E.1 (Main Study, Unrestricted)

tab_modelFit_s2 %>%
  filter(variable != "(Intercept)") %>%
  mutate(
    Variable = recode_factor(variable, disadv_nuke = "eta Nuclear", adv_chance = "alpha1 Chance of success", adv_mil = "alpha2 Prevents casualties", 
                             disadv_ally = "delta1 Allies disapprove",
                             disadv_civ = "delta2 Civilian casualties", disadv_env = "delta3 Environmental damage"),
    Estimate = paste0(format(round(estimate, 2)), " (", format(round(se, 2)), ")")
  ) %>%
  arrange(Variable) %>%
  select(Variable, Estimate) %>%
  xtable() %>%
  print(only.contents=T, hline.after = nrow(.),include.rownames=F, include.colnames=F) %>%
  gsub("eta", "$\\\\eta$~ & ", .)%>%
  gsub("alpha1", "$\\\\beta_1$~~ & ", .)%>%
  gsub("alpha2", "$\\\\beta_2$~~ & ", .)%>%
  #  gsub("alpha3", "$\\\\beta_3$~~ & ", .)%>%
  gsub("delta1", "$\\\\beta_3$~~ & ", .)%>%
  gsub("delta2", "$\\\\beta_4$~~ & ", .)%>%
  gsub("delta3", "$\\\\beta_5$~~ & ", .)%>%
  write("tab_model_s2.txt")

#Panel C of Table E.1 ("Pilot, restricted model")
tab_modelFit_restrict_s1 %>%
  filter(variable != "(Intercept)") %>%
  mutate(
    Variable = recode_factor(variable, disadv_nuke = "eta Nuclear", ovr_adv = "blah Net advantages"),
    Estimate = paste0(format(round(estimate, 2)), " (", format(round(se, 2)), ")")
  ) %>%
  arrange(Variable) %>%
  select(Variable, Estimate) %>%
  xtable() %>%
  print(only.contents=T, hline.after = nrow(.),include.rownames=F, include.colnames=F) %>%
  gsub("eta", "$\\\\eta$~ & ", .)%>%
  gsub("blah", "$\\\\beta$~~ & ", .)%>%
  write("tab_model_s1_restrict.txt")

#For Panel d (Main study, restricted model)
tab_modelFit_restrict_s2 %>%
  filter(variable != "(Intercept)") %>%
  mutate(
    Variable = recode_factor(variable, disadv_nuke = "eta Nuclear", ovr_adv = "blah Net advantages"),
    Estimate = paste0(format(round(estimate, 2)), " (", format(round(se, 2)), ")")
  ) %>%
  arrange(Variable) %>%
  select(Variable, Estimate) %>%
  xtable() %>%
  print(only.contents=T, hline.after = nrow(.),include.rownames=F, include.colnames=F) %>%
  gsub("eta", "$\\\\eta$~ & ", .)%>%
  gsub("blah", "$\\\\beta$~~ & ", .)%>%
  write("tab_model_s2_restrict.txt")


#   APPENDIX: Data fit the model ------------------
#This will make Figure E.1

tab_pred = tab_plot_ovr_index %>%
  group_by() %>%
  select(Study, types) %>%
  full_join(expand.grid(Study = c("Study 1", "Study 2"),ovr_adv = ((-500:500)/100))) %>%
  filter(!(abs(ovr_adv)>4 & Study=="Study 1")) %>%
  mutate(
    latent = case_when(
      Study == "Study 1" ~ tab_modelFit_restrict_s1$estimate[3]*ovr_adv + tab_modelFit_restrict_s1$estimate[2]*grepl("Nuclear", types),
      Study == "Study 2" ~ tab_modelFit_restrict_s2$estimate[3]*ovr_adv + tab_modelFit_restrict_s1$estimate[2]*grepl("Nuclear", types)
    ),
    estimate = invlogit(latent),
    conf.low = NA_real_,
    conf.high = NA_real_,
    Study = recode(Study, `Study 1`="Pilot", `Study 2`="Main Study")
  ) %>%
  arrange(ovr_adv)

g2 = g + geom_path(data = tab_pred%>%
                     filter(types %in% c("Nuclear vs. Conventional", "Conventional vs. Conventional")), alpha=.35, size=2)

ggsave("pappendix_modelPred_conjoint.pdf", g2, "pdf", width=6.5, height=3.5)
ggsave("FigE.1.tiff", g2, "tiff", width=6.5, height=3.5)

log_close()
